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ABSTRACT 


This paper contains the method used to solve Incompressible Navier Stokes 
(INS) equation. This study is conducted to see how the methods used in this 
paper gives a solution over a particular domain. Fraction Step Method is 
used for solving INS equation. Non-dimensional equation is taken for 
formulation. It is solved in two parts. With first part, intermediate velocity 
is calculated. Then next part is used to calculate pressure. Once pressure is 
obtained, then velocity is corrected for next time step. This is called 
Pressure correction method. This process is continued until steady state 
solution is obtained. After that, vorticity is calculated, and plotted. These 
plots are analyzed and compared with the experimental solution. Further a 
conclusion is made regarding whole paper. 


INTRODUCTION 

Property of a fluid element can be predicted for a domain 
by using the following three governing equation: 

> Continuity Equation: Based on conservation of mass. 

> Momentum Equation: Based on Newton's second law 
of motion. 

> Energy Equation: Based on conservation of energy 
within a system. 

These three governing equations were very accurate in 
determining the flow properties under a particular 
circumstance. Depending upon conditions, these equations 
are formulated for a specific problem. Like solving for 
compressible, incompressible, steady, unsteady, laminar, 
turbulent, isothermal flows. All the governing equations 
are partial differential equations therefore solved using 
numerical technique using power of computer. 

Incompressible Navier Stokes equation is a case where 
spatial derivative for density is equal to zero. Thus, it gives 
a very simple equation for formulation. Non-dimensional 
form of Incompressible Navier Stokes equation will be 
used for the formulation: 

^ + (u*. V‘)u* = -w* + 

Here first term on LHS is the local derivative of velocity, 
second term is convective term for velocity, on RHS, 
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gradient of pressure is taken first and Laplacian of 
velocity, also called diffusion term. There is one more term 
which is ignored is gravitational acceleration term. 

There is one approach for formulating Incompressible 
Navier Stokes Equation by formulating a general case for 
Compressible flows, and then put all spatial density 
derivatives as zero. But there is an error in this approach. 
For compressible flows, time step is calculated by the 
given formula: 

1 

At < - - 

\u\ , bl , _ 1 i ! ~ 

Ax^ Ay (Ayy 

For compressible flows, 'a' (speed of sound) is finite. So, it 
will give a finite value of time step. But for incompressible 
flows, the speed of sound is theoretically infinite. Which 
will give time step At = 0. This is not desirable as it not 
even start the calculation. It is also seen as when a 
compressible CFD solver is run for a Mach number near to 
zero, it takes more time steps to converge the solution. 
And most of the time it has tendency to be unstable at such 
low Mach numbers. Therefore, a different approach is 
used to make incompressible solver. Finite difference 
method is used for generating numerical solution of 
Incompressible Navier Stokes equation. FORTRAN 
Language will be used for the formulation of this solver. 
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Figure 1: Grid used for formulation of Incompressible 
Navier Stokes Equation. 

Above figure contains a rectangular grid which represents 
the domain for this solver. All the small rectangular boxes 
are representing fluid element. Property of fluid element is 
calculated at the intersection point of two elements, also 
called as nodes. Length of a fluid element in X direction is 
Ax while the same in Y direction is Ay. 


Fraction step method: 

For solving the non-dimensional incompressible equation, 
fraction step method is used. For compressible flows 
pressure can be calculated using equation of state. But in 
incompressible form, equation of state cannot be used. 

^+(u.V)u^-Vp + ^J^u ( 1 ) 


For solving Incompressible Navier Stokes equation 
pressure gradient is needed. So, above equation is solved 
in two parts: 


= ( 2 ) 


du 


= -Vp 


( 3 ) 


Motive is to get velocity of fluid at next time step, so first 
velocity is calculated without taking the effect of pressure 
gradient. This velocity obtained is not the exact velocity 
for the fluid element at next time step. Therefore, this 
obtained velocity is called intermediate velocity 'u'. This 
intermediate velocity is used in next equation to calculate 
pressure. Once pressure is calculated then with the help of 
pressure, velocity is corrected and new velocity is 
calculated for next time step. 


Calculation of Intermediate Velocity(u): 

For calculating intermediate velocity, eqn. 2 is used. 


^+(u.V)u^—V^u 

dt ^ ^ Re 


First term in LHS is the local derivative of velocity, next is 
convective term for velocity. The term on RHS is inverse of 
Reynolds Number multiplied by Laplacian of velocity 
vector. 


here 'n' is the current time step. 

Fort the convective term, upwind method is used. This is 
so because accuracy is not the priority but calculation 
time. Upwind method is therefore preferred for 
calculation of convection term. 

Upwind scheme is: 

(or) 

If u which is outside of the bracket is positive then 
backward difference is used, if it is negative then forward 
difference is used. 

For diffusive term on RHS, which is Laplacian of velocity 
vector, second order central difference is used: 

t ^^2 J t ^y2 J 

After arranging all the terms, intermediate velocity is 
calculated. 


Calculation of Pressure: 

Pressure for fluid elements is calculated using equation 3: 
du ^ 


First order finite difference is used for formulating local 
velocity derivative. For pressure term, central difference is 
used, because it converges and chances of instability is 
very least compared to forward and backward difference: 

(n+l) _ 

At ^ 2AX ' 


Above equation cannot be solved because there are two 
unknowns, one is and second is 'p'. This can be 

solved with property of incompressible flows. It states that 
divergence of velocity is zero. If divergence of above 
equation is taken then will be reduced to zero. 

Intermediate velocity cannot be zero because it was not 
calculated taking the effect of pressure: 


This gives a Pressure Poisson equation: 



Above equation contains Laplacian of pressure on LHS, 
which is equal to the source term on RHS. This source 
term is divergence of intermediate velocity over time step. 
This Pressure Poisson Equation is then solved in the given 
domain to obtain the value of pressure. 


Solving Pressure Poisson Equation: 

For solving Pressure Poisson Equation, it is converted into 
linear algebraic equation. For Laplacian of pressure, 
second order central difference is used: 


72^ _ r ^t+1 J')^(t-1 J') ij) _L r ^(ij+i)+P(ij-i) 






^} + {- 


Ay^ 


Local velocity derivative is calculated with first order 
forward difference method: 
du _ u-u^'^^ 
dt At 


For source term, first order forward difference is used: 

^ AX ^ Ay ^ 
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Figure 2: Grid including ghost points for Boundary 
condition of Pressure Poisson Equation 

For starting the calculation boundary conditions are 
needed. Here Neumann boundary condition is used which 
states that normal component of pressure is equal to 
constant. Here the constant terms take is zero. 


So: ^ = 0 for x = 0, n; and ^ = 0 for y = 0, m 

ox oy 


First order central difference is used to formulate above 

conditions: 

P(i+ij)-P(i-ij) _ 


2 AX 


0 for i = 1, n 


And 

P(ij+l)“P(ij-l) ci c ■ 'I 

^^ = 0 for 1 = 1, m 
2Ay ^ 


Above conditions are used with ghost points and then 
substituted in Pressure Poisson Equation for each 
boundary condition. Ghost points are at j = 0 , m+1 and i = 
0, n+1. 

Above finite equations are simplified for boundary points: 
Fori = 1 ^p[ 2 ,j) = p[o,j); 

For i = n ^ P[n-i,D = P[n+i,j) 

For j = 1 ^ P(i, 2 )=p[i,o) 

For j = m ^ P(i,m-l)=P[i,m+l) 

Substituting these conditions in Pressure Poisson 
Equation, solution at boundary points are calculated by 
following equations: 


P(ij) 




2 (- 


for i = 1 




P(n,;) “ 


L J L J J 


2 (—2 - 2 ) 


for i = n 


p(u) =-^ 


Vii.ra) 


^ Ax^ ^ Ay' 


^ ■' ^At ■' 


2 (- 


for j = m 


Ax^ Ay^ 


After getting pressure on boundary points, Pressure Poisson equation is solved on internal points: 


{■ 


P(i+l,j) + P(i-l,j)~, , rP(U+l) + P(iJ-l) 


Pa.n 


A 


^} + {- 


A 




^x(i+ij) . ,r^y(i ,7 + 1 ) ^y(ij) 


A X 


^) + (- 


A y 


■)} 


2 (- 


A X^ 


■ + 


A y 


r) 


Using Gauss-Seidel method, whole Pressure Poisson Equation is iterated until solution converges. At last a set of values 
will be obtained which will be used to calculate velocity at next time step. 


Calculation of Velocity at next time step: 

After getting pressure at all the points, the effect 
pressure is used to calculate velocity for next time step: 

.gy) = -A 


of 


ForX: 
ForY: u, 


(y) 


2Ay 


This method for using the pressure to correct and obtain 
next time step velocity is also called Pressure Correction 
Method. To obtain steady state solution, continuous 
iterations are done using obtained values until the 
solution converges for velocity. 

Result and Analysis: 

After the solver is converged, a set of data is obtained from 
the solution. This data contains Pressure, X velocity and Y 
velocity. There is an additional step in the program, for 
analyzing the results, vorticity is calculated. This is done 
by: 

Vorticity ■ 


dv 

dx 


du 

dy 


Above equation is formulated using first order finite 
difference method. The choice between forward and 
backward scheme is made depending on the boundary 


points. For all the points, forward difference is used but at 
i = n and j = m, backward difference is used. 

For test case, condition of lid-driven cavity is simulated. 
Boundary condition for velocity is zero on all boundaries 
except at top boundary, where X velocity is taken as one. 
Test case is run for Reynolds Number equal to 100, 500 
and 1000. A square grid of 100 by 100 is taken for test 
case. Results are obtained with a time step of 0.004 for all 
three conditions. Code was run was time period of 60 secs 
to get very accurate results and ensure steady state 
solution is obtained. This particular time step is chosen 
because it was the maximum value at which solution was 
stable. If time step is increased any further then there was 
instability in the code. 

Results show that there is a significant variance from 
Reynolds Number 100 to 500. After that that change is not 
that prominent. For Reynolds Number 100: Formation of a 
thick boundary layer is seen across the side walls. The 
velocity at the bottom of the contour decreases 
significantly. Despite this, there are formation of strong 
vortices around the top corners of the domain. Vortices 
pass through the center for lower Reynolds Number. For 
Reynolds Number 500 to 1000: A significant change is 
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observed from Re = 100. Velocity is more around the 
center and top boundary of the domain compare to other 
regions. Unlike at Re 100, where velocity gradually 
decreases as going to bottom boundary. It is also observed 
that as Re No. increases, a thin boundary layer forms at the 
side walls. Strong vortices are formed at the top corners, 
and here vortices do not pass through the center. Because 
velocity is more at the center, flow is circulating about the 
center of the domain. As Re. No. increases, this circulation 
moves further away from the center. For Re 1000, 
boundary layer is thinner, and there is formation of small 
vortices at the bottom. 


Vorticity at Re 100 
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Figure 3: Vorticity Contour Plot for Reynolds no: 100 



Vorticity at Re 500 
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Figure 4: Vorticity Contour Plot for Reynolds no: 500 
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Figure 4: Vorticity Contour Plot for Reynolds no: 1000 


Vorticity at Re 1000 
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Center of the primary eddy remains constant below the 
top lid for Reynolds number 100. While for higher Re. No. 
it lies at the center of the domain. One common 
observation for all the cases is boundary layer is thinner at 
top boundary because of X component of velocity there. 

Conclusion: 

The solver for Incompressible Navier Stokes equation was 
successfully formulated. Use of finite difference was made 
and how different schemes affect the solution of the 
problem was observed. A small change in input 
parameters can result to instability in the solution and 
output inaccurate or no result. Formulation of Pressure 
Poisson Equation and boundary condition were seen as 
the most critical part of the solver. A very simple error of 
array dimension misplacing can diverge the solution. For 
different set of numerical schemes, different solutions are 
obtained. If forward or backward difference is used in 
pressure term while correcting velocity for next time step, 
then solution never converges and causes instability. This 
shows that proper set of schemes are used to give stable 
results. After obtaining solution, data were plotted in 
OriginLAb pro software for visualization and analysis. It 
was seen how the dynamics of flow differs as Reynolds 
Number increases. Further then the obtained solution was 
compared with standard data available on internet and 
accuracy was obtained while obtaining the solution. 
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